Lattice susceptibility for 2D Hubbard Model within dual fermion method 
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In this paper, we present details of the dual fermion (DF) method to study the non-local correction 
to single site DMFT. The DMFT two-particle Green's function is calculated using continuous time 
quantum monte carlo (CT-QMC) method. The momentum dependence of the vertex function is 
analyzed and its renormalization based on the Bethe-Salpeter equation is performed in particle-hole 
channel. We found a magnetic instability in both the dual and the lattice fermions. The lattice 
fermion susceptibility is calculated at finite temperature in this method and also in another recently 
proposed method, namely dynamical vertex approximation (DFA). The comparison between these 
two methods are presented in both weak and strong coupling region. Compared to the susceptibility 
from quantum monte carlo (QMC) simulation, both of them gave satisfied results. 

PACS numbers: 71.10.Fd 



I. INTRODUCTION 

Strongly correlated electron systems, such as heavy 
fermion compounds, high-temperature superconductors, 
have gained much attention from both theoretical and 
experimental point of view. The competition between 
the kinetic energy and strong Coulomb interaction of 
fermions generates a lot of fascinating phenomena. Vari- 
ous theoretical approaches have been developed to treat 
the regime of intermediate coupling. The widely used 
perturbative methods, such as random phase approxi- 
mation (RPA), fluctuation exchange (FLEX)ii^, and the 
two-particle self-consistent (TPSC)^i^ method are based 
on the expansion in the Coulomb interaction which is 
only valid in weak-coupling. To go beyond the perturba- 
tive approximation and to gain insight of the correlation 
effects of the fermion systems, new theoretical methods 
are needed. Dynamical mean field theory (DMFT)^'^'^ is 
a big step forward in the understanding Metal-Insulator 
transition. 

Dynamical mean field theory maps a many-body inter- 
acting system on a lattice onto a single impurity embed- 
ded in a non-interacting bath. Such a mapping becomes 
exact in the limit of infinite coordination number. All 
local temporal fluctuations are taken into account in this 
theory, however spatial fluctuations are treated on the 
mean field level. DMFT has been proven a successful 
theory describing the basic physics of the Mott-Hubbard 
transition. But the non-local correlation effect can't al- 
ways be omitted. Although, straight forward extensions 
of DMFT— i^iiSiiiii^ have captured the influence of short- 
range correlation, these methods are still not capable of 
describing the collective behavior, e.g. spin wave excita- 
tions of many-body system. At the same time, most of 
the numerically exact impurity solvers require a substan- 
tial amount of time to achieve a desired accuracy even on 
a small cluster, which makes the investigation of larger 
lattice to be impossible. 

Recently, some efforts have been made to take 
the spatial fluctuations into account in different 
ways^^i^i^'^i^. All these methods construct the non- 



local contribution of DMFT from the local two-particle 
vertex. The electron self-energy is expressed as a function 
of the two-particle vertex and the single-particle propa- 
gator. The cluster extention of DMFT considers the cor- 
relation within the small cluster. Compared to these, the 
diagrammatic re-summation technique involved in these 
new methods makes them only approximately include the 
non-local corrections. While, long range correlations are 
also considered in these methods and the computational 
burden is not serious. 

In this paper we will apply the method of Rubtso"\J^ 
to consider the vertex renormalization of the DF through 
the Bethe-Salpeter equation. Lattice susceptibility is cal- 
culated from the renormalized DF vertex. 

The paper is organized as follows: In Sec. |TT]we sum- 
marize the basic idea of the DF method and give details 
of the calculation. The DMFT two-particle Green's func- 
tion and the corresponding vertex calculation are imple- 
mented in CT-QMC in Sec. IIIII The frequency depen- 
dent vertex is modified through the Bethe-Salpeter Equa- 
tion to obtain the momentum dependence in Sec. IIVI In 
Sec. |V] we present the calculation of the lattice suscep- 
tibility and compare it with QMC results and also the 
works from Toschii^. The conclusions are summarized in 
Sec. IVIl where we also present possible application. 



II. THE DF METHOD 

We study the general one-band Hubbard model at two 
dimensions 



^ = H efc.acLcfc,, + UY^ n,^n,^ 
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Cktr) creates (annihilates) an electron with spin-tr 
momentum k. The dispersion relation is e^ = 
—'2tJ2i=i cos ki, where N is the number of lattice sites. 
The basic idea of the DF method^"' is to transform the 
hopping between different sites into coupling to an aux- 
iliary field /(/^). By doing so, each lattice site can be 



viewed as an isolated impurity. The interacting lattice 
problem is reduced to solving a multi-impurity problem 
which couples to the auxiliary field. This can be done 
using the standard DMFT calculation. After integrating 
out the lattice fermions c{c^) one can obtain an effective 
theory of the auxiliary field where DMFT serves as an 
starting point of the expansion over the coupling between 
each impurity site with the auxiliary field. 

To explicitly demonstrate the above idea we start from 
the action of DMFT which can be written as 



5[c+,, 
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the DMFT self consistency. Therefore the first non-local 
contribution is given by diagram (b). The self-energy for 
these two diagrams are 

S«(fci) - -§ E Gi,{k,Wf^, {!.,,.' ;,.',:.) (5a) 
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where A^, is the hybridization function of the impurity 
problem defined by Sl^^ which is the action of an iso- 
lated impurity at site i with the local Green's function 
gu- Using the Gaussian identity, we decouple the lattice 
sites into many impurities which couple only to the field 
/ 

S[c\c;P,f] = E^™P+E[3-"'(4./fe.. + ^-c.) 

i k.u,a 

+9-\A. - efe)- VL4.J (3) 

The equivalence of Eqs. ([2]) and Q form an exact rela- 
tion between the Green's funtion of the lattice electrons 
and the DF. 



G,.fe = g-''{A, - Ek) 






(A^ - cfc) 



(4) 



This relation is easily derived by considering the deriva- 
tive over efc in the two actions. Eq. (U) allows now to 
solve the many-body "lattice" problem based on DMFT 
which is different from the straight forward cluster exten- 
sion. The problem is now to solve the Green's function 
of the DF Gf, ). . It is determined by integrating Eq. ([3]) 
over c' and c yielding a Taylor expansion series in powers 
of /^ and /. The Grassmann integral ensures that / and 
/ appear only in pairs associated with the lattice fermion 
n-particle vertex obtained from the single-site DMFT cal- 
culation. In this paper we restrict our considerations to 
the two-particle vertex j^'^' . 

Expanding the Luttinger-Ward functional in 7^, the 
first two contributions to the self energy function are the 
diagrams shown in Fig. [TJ Diagram (a) vanishes for 
the bare DF since this diagram exactly corresponds to 
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FIG. 1: The first two self-energy diagrams. They are com- 
posed of the local vertices function and DF propagator. 



Here space-time notation is used, k — [k^v), q — {q,Lu). 
Fermionic Matsubara frequency is Vn — (2n -t- 1)tt/P, 
bosonic frequency is LUm = 2rm: / (3 where [3 is the inverse 
temperature. Together with the bare DF Green's func- 
tion G'o(/c) — —g1l\{Ay - Cfe)"^ -I- g^], the new Green's 
function can be derived from the Dyson equation 

[G''{k)]-' = [G^,{k)]-'-^\k) (6) 

The algorithm of the whole calculation is: 

1. Set initial value of A^ for the first DMFT loop. 

2. Determine the single-site DMFT Green's function 
gi, from the hybridization function A^. The self- 
consistency condition ensures that the first diagram 
of the DF self-energy is very small. 

3. Go through the DMFT loop once again to calcu- 
late the two-particle Green's function and corre- 
sponding 7-function. The method for determining 
the 7-function is implemented for both strong and 
weak-coupling CT-QMC in the next section of this 
paper. 

4. Start an inner loop calculation to determine the DF 
Green's function and in the end the lattice Green's 
function. 

(a) From Eqs. ([5a)) . (|5bll and the Dyson equation 
dHI) to calculate the self-energy of the DF. 

(b) Repeatly use Eq. jSa]), (l5b| and Eq. jB]) until 
the convergence of the DF Green's function is 
achived. 

(c) The lattice Green's function is then given by 
Eq. (gl) from that of the DF. 

5. Fourier transform the momentum lattice Green's 
function into real space. And from the on-site com- 
ponent Gii to determine a new hybridization func- 
tion Aj, which is given by Eq. ^. 

6. Go back to the Step 3. and iteratively perform 
the outer loop until the hybridization A^, doesn't 
change any more. 



Although diagram (a) is exactly zero for the bare DF 
Green's function, it gives non-zero contribution from the 
second loop where the DF Green's function is updated 
from Eq. ^. As a result, the hybridization function 
should also be updated before the next DMFT loop is 
performed . This is simply done by setting the local full 
DF Green's function to zero, together with the condition 
that the old hybrization function forces the bare local DF 
Green's function to be zero {J^k ^Jk = 0)i we obtain a 
set of equations 



1 Y,[G.,k - (Af-' - 6,)-%^(Ar'" - ^kf = (7a) 
fc 

\ Y,[Gi^, - (Ar - ^.)-%^(Ar - ^k? = (7b) 
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which yields 
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This equation finally gives us the relation between the 
new and old hybridization function. 
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(9) 



In the whole calculation, the DF perturbation calcula- 
tion converges quickly. The most time consuming part of 
this method is the DMFT calculation of the two particle 
Green's function. There are some useful symmetries to 
accelerate the calculation. As already pointed ou t^^'^^ , 
it is convenient to take the symmetric form of the inter- 
action term. The two particle Green's function is then a 
fully antisymmetric function. Such fully antisymmetric 
form is very useful to speed up the calculation of the two 
particle Green's function. One does not need to calculate 
all the frequency points within the cutoff in Mastsubara 
space, a few special points are calculated and the values 
for the other points are given by that of those special 
points through antisymmetric property. In the DF self 
energy calculation, we always have the convolution type 
of momentum summation which is very easy to be calcu- 
lated by fast fourier transform (FFT). 



III. CT-QMC AND TWO-PARTICLE VERTEX 

From the above analysis, the key idea of the DF 
method is to construct the nonlocal contribution from the 
auxiliary field and the DMFT two-particle Green's func- 
tion. Therefore it is quite important to accurately deter- 
mine the two-particle vertex. Here we adapt the newly 
developed CT-QMC method^^i^iiSa to calculate the two 
particle Green's function x- 

First we briefly outline the CT-QMC technique. For 
more details, we refer the readers to ^°i^^'^^ . Here we dis- 
cuss the two-particle Green's function and some numeri- 
cal implemetations in more detailed. Two variants of the 



CT-QMC methods have been proposed based on the di- 
agrammatic expansion. Unlike the Hirsch-Fye method, 
these methods don't have a Trotter error and can ap- 
proach the low temperature region easily. In the weak- 
coupling method^° the non-interacting part of the parti- 
tion function is kept and expanded the interaction term 
into Taylor series. Wick's theorem ensures that the cor- 
responding expansion can be written into a determinant 
at each order 



E 






dri ■ • ■ drfcC"*" dei[D^Di] (10) 



with 



D^D^^ 



Gt(ri -Tfe) 



Gi(Tfe-ri) 



(11) 
where 5*0 is the non-interacting action and G° is the Weiss 
field, and the one-particle Green's function is measused 
as 

G{v) = G\v) ~ iG°(z.)^Af,,e-(---^)G"(z.) (12) 



/3 



*j 



In the strong coupling method the effective action is ex- 
panded in the hybridization function by integrating over 
the non-interacting bath degrees of freedom. Such an 
expansion also yields a determinant. 

Z ^ TrTre-^'- H E ^ / '^^i ' ' ' ^^t / ^rf ■ ■ ■ dr^ 

a- ka 

I A(Tf-rf) ... A(Tf-r,!j \ 
^.{r^)\ ... ■•. ... *i(r^)(13) 

Here *(t) — (ci(t), C2(t), . . . ,Ck^(T))- The action is 
evaluated by a Monte Carlo random walk in the space 
of expansion order k. Therefore the corresponding hy- 
bridization matrix changes in every Monte Carlo step. 
One particle Green's function is measured from the ex- 
pansion of hybridization function as G{t!^ — r*) = Af^.j. 
M is the inverse matrix of the hybridization function. 
Apparently one needs to calculate this inverse matrix in 
every update step which is time consuming, fortunately 
it can be obtained by the fast-update algorithm^S. 

At the same time such a relation allows direct mea- 
surement of the Matsubara Green's function 



G{i 



^^' 



•M 



«j^ 



(14) 



Compared with the imaginary time measurement, it 
seems additional computational time is needed for the 
sum over every matrix elements Mij. K. Haule proposed 
to implement such measurement in every fast update pro- 
cedure which makes sure that only linear amount of time 
is needed22i. 



In our calculation the Green's function is measured 
in the weak-coupling CT-QMC at each accepted up- 
date which greatly reduces the computational time. The 
weak-coupling CT-QMC normally yields a higher per- 
turbation order k than the strong-coupling CT-QMC. It 
seems that the performance of the strong-coupling CT- 
QMC is better^^. Concerning the convergence speed, the 
weak-coupling CT-QMC is almost same as the strong- 
coupling one under the above implementation together 
with a proper choice of a, since in strong-coupling CT- 
QMC more Monte Carlo steps are needed usually in order 
to smooth the noise of Green's function at imaginary time 
around /3/2 or at large Matsubara frequency points. Fur- 
thermore, the weak-coupling CT-QMC is much easier im- 
plemented for large cluster DMFT calculation, in which 
case the strong-coupling method needs to handle a big 
eigenspace. In this paper we mainly use weak-coupling 
CT-QMC as impurity solver, while all the results can 
be obtained in the strong-coupling CT-QMC which was 
used as an accuracy check. 

Similarly, we adapt K. Haule's implementation to cal- 
culate the two-particle Green's function in frequency 
space. In the weak coupling CT-QMC, the non- 
interacting action has Gaussian form which ensures the 
applicability of Wick's theorem for measuring the two 
particle Green's function 



Xcrcr'{vi,1^2,1^3,l'i) = T[Ga{vi, V2)Ga' {vs^Vi) 



- 5aa,Ga{vi,Vi)Ga{v3„V2)\{l^) 

The over-line indicates the Monte Carlo average. In 
each Monte Carlo measurement, G{v, v') depends on two 
different argument v and u' , only in the average level, 
G(y, v') ~ G{v)5^^jji is a function of single frequency. 
In each fast-update procedure, the new and old G(y, v') 
have a closed relation which ensures that one can deter- 
mine the updated Green's function G^'^'^{v, v') from the 
old one G'-'^'''(y,v'). For example, adding pair of kinks 
and supposing before updating the perturbation order is 
k, then it is k-\-l for the new M-matrix. The new inserted 
pair is at fc -|- 1 row and fc -I- 1 column. 

G^'''"{v,v')-G°^'^{v,v') 



j^jNew 



j^^^G"{iy) {xL XR-XR- 6-'"^+^ 



/3 
-XL ■ e"''^^^+^ + e-'''^^-+i+'"''^^-+i| G'"(i/'X16) 

Here, XL = ELi e"""^*'-^*, XR = Y!]=i e'""' ''^ Rj and 
Li, Rj have the same definition as in Re^. In every step, 
one only needs to calculates the Green's function when 
the update is accepted and only a few calculations are 
needed. A similar procedure for removing pairs, shift- 
ting end-point operation can be used. Such method is 
also applicable in the segment picture of strong-coupling 
CT-QMC. In the weak-couphng CT-QMC, such an im- 
plementation greatly improves the calculating speed in 
low temperature and strong interaction regime^*^. Once 



one obtains the two frequency dependent Green's func- 
tion in every monte carlo step, the two-particle Green's 
function can be determined easily from Eq. (|15p . The 
two-particle vertex is then given from the following equa- 
tion: 

,-(..') =,,5¥^^^;^^#^^ (17) 



ga{v)ga{v + w)ff<T'(^' + '^)9<y'{v') 



where 



(18) 
is the bare susceptibility. For the multi-particle Green's 
function, it still can be constructed from the two fre- 
quency dependent Green's function G{v,v')^ but more 
terms appear from Wicks theorem. Simply, when set 
V = v' one can calculate the one-particle Green's funtion 
easily. 



IV. MOMENTUM DEPENDECE OF VERTEX 

As mentioned earlier diagram (a) in Fig. [l]only gives 
the local contribution. The first non-local correction in 
the DF method is from diagram (b). Momentum de- 
pendences comes into this theory through the bubble- 
like diagram between the two vertices which yields the 
momentum dependence of the DF vertex. The natural 
way to renormalize vertex is through the Bethe-Salpeter 
equation. Since the DMFT vertex is only a function of 
Matsubara frequency, the integral over internal momen- 
tum k and k' ensures that the full vertex only depends 
on the center of mass momentum Q. The Bethe-Salpeter 
equation in the particle-hole channe l^^'^^ are shown in 
Fig. H 

From the construction of the DF method, we know 
the interaction of the DF is coming from the two parti- 
cle vertex of lattice fermion which is obtained through 
DMFT calculation. In the Bethe-Salpeter equation, it 
plays the role of the building-block. The corresponding 
Bethe-Salpeter equation for these two channels are 



T 



^ E l7" {yy)G\k")G\k" + Q)V^^''^""'^' {v\v') 

(19a) 
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N 



E 7r(^, ,^")G'ik")G''ik" + Q)r^''i''^'^(^", z.') 



(19b) 



Here, the short hand notation of spin configuration is 
used. 7'^'^ represents 7°'°'°' °" , while 7'^'^'^°' is denoted by 
7°'°' where a = —a. Tp'^^^p'^^'' are the full vertices in the 
Sz — and S'z = ±1 channel, respectively. G'^ is the full 
DF Green's function obtained from section [H] which is 







a 
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FIG. 2: S^ = (phO) and S^ = ±1 (phi) particle-hole chan- 
nels of the DF vertex, between vertices there are two full 
DF Green's function. The Sz = ±1 component is the triplet 
channel, while that for Sz — can be either singlet or triplet. 



kept unchanged in the calculation of the Bethe-Salpeter 
Equation. Different from the work of S. Brener— , we 
solve the above equations directly in momentum space 
with the advantage that in this way we can calculate the 
susceptibility for any specific center of mass momentum 
Q and it's convenient to use FFT for investigating larger 
lattice. 

In Eq. p9|) one has to sum over the internal spin 
indices in the Sz = channel which is not present in 
5'z = ±1 channel. One can decouple the Sz = chan- 
nel into the charge and spin channels 'yc{s) — 1°'°' i 7°'°' 
which can be solved seperately, and it turns out that the 
spin channel vertex function is exactly same as the that 
in Sz = ±1 channel, see e.g. P. Nozieres"'^^. Such relation 
is true for the DMFT vertex, and was also verified for 
the momentum dependent vertex in the DF method^^. 
In our calculation, we have solved the Sz — channel by 
decoupling it to the charge and spin channel, while the 
phi channel is not used. 

Once the converged momentum dependent DF vertex 
is obtained, one can determine the corresponding DF sus- 
ceptibility in the standard way by attaching four Green's 
functions to the DF vertex. 



xr'(Q) = x"(Q)^ 
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(20a) 



xf(Q) = xS(Q)+ 

|J ^ Gi{k)Gi{k + Q)r'^''{Q)Gi{k')Gi{k' + Q) 



iV2 



k.k' 



(20b) 



The momentum sum over k and k' can be performed 
independently by FFT becasue the DF vertx V^"^' {Q) 
only depends on the center of mass momentum Q. 

Now the z-component DF spin susceptibility (S'^-5^) = 
hiXd ~ Xdl ^^^ ^^ determined from the spin channel 
component calculated above. In Fig. [31 x^^ = x^^ — Xo^ 
is shown for U/t = 4 at temperatures (3t — 4.0 (left 
panel) and f3t = 1.0 (right panel). With the lowing down 
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FIG. 3: The nontrivial part of the DF spin susceptibilities as 
a function of momentum in 2D Hubbard Model at U/t = 4.0, 
Pt = 1.0 (right panel) and f3t = 4.0 (left panel). Here 32 x 
32 momentum points are used in the first Brillouin zone. 



of temperature the DF susceptibility grows up, especially 
at wave vector (tt, tt). The momentum k^ and ky run from 
to 2tt. The susceptibility is strongly peaked at the wave 
vector (tt, tt) at the low temperature case and the peak 
value becomes higher and higher. The magnetic insta- 
bility of the DF system is indicated by the enhancement 
of the DF susceptiblity. The effect of momentum depen- 
dence of vertex is clearly visible in this diagram. The bare 
vertex which is only a function of frequency becomes mo- 
mentum dependent through the Bethe-Salpeter equation. 
Later on we will see that such momentum dependent ver- 
tex plays a very important role in the calculation of the 
lattice fermion susceptibility. 



V. LATTICE SUSCEPTIBILITY 

The strong antiferromagnetic fluctuation in 2D system 
is indicated by the enhancement of the DF susceptibil- 
ity at the wave vector (tt, tt) shown in Fig [31 This is 
the consequence of the deep relation between the the 
Green's function of the lattice and the DF, see Eq. Q. 
In order to observe the magnetic instability of the lat- 
tice fermion directly, we have calculated the lattice sus- 
ceptibility based on the DF method. By differentiating 
the partition function in Eqns. ^ [3]) twice over the ki- 
netic energy, we obtain an exact relation between the 
susceptibility of DF and lattice fcrmions. After some 
simplifications^^, it is given by 

Xf{Q)^x'f{Q) + 

^'. J2 G'{k)G'{k + Q)r$(i., v')G'{k')G'{k' + qn) 



N^ 



k,k' 



Here G" cannt be interpreted as a particle propagator, it 
is defined as: 



G'(fc) 



G'^(fc) 



g4A,-e{k)] 



(22) 



Again, the sum is performed over internal momentum 
and frequency fc, k' which is performed by FFT and rough 
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FIG. 4: The uniform spin suscetibility of the DF using the 
bare vertex (only frequency dependent) and the full ver- 
tex(vertex from the Bethe-Salpeter quqation) for half filled 
2D Hubbard model at U /t = 4.0 and various temperatures. 
These results reproduce the similiar solution in comparison 
with the calculation of finite size of QMC. 



sumniing over a few Matsubara points. Again as in Eq. 
(|4|), this equation established a connection between the 
lattice susceptibility and the DF susceptibility. From this 
point of view, it is easy to understand that the instability 
of DFs generates the instability of the lattice fermions. 

One can also find relations for the higher order Green's 
function of the DF and the lattice fermions in the same 
way. This emphasizes the similar nature of the DF and 
lattice fermions except that DF possess only non-local in- 
formation, since the DMFT self-consistency ensures that 
the local DF Green's function is exactly zero. 

The lattice magnetic susceptibility is calculated using 
the following definition 

= 2{xy-xy) (23) 

where XfihT) = ([wj^tW "'^vi W] ^ ["-o,t(0) -rioa(O)])- 
Xf represents the lattice susceptibility in order to distin- 
guish with that of the DF. 

We have used two different ways to calculate the lattice 
susceptibility. First we have solved the above equation 
using the bare vertex T(i>, J^'; w) which is obtained from 
the DMFT calculation. In contrast, the second calcula- 
tion was performed using the full DF vertex. In both 
of these calculations, the full one particle DF Green's 
function was used. The momentum dependent of the DF 
vertex is obtained through the calculation of the Bethc- 
Salpeter equation. The lattice susceptibility is expected 
to be improved if we use the momentum dependence DF 
vertex. In this way, we can understand the effect of mo- 
mentum dependence in the DF vertex. 

In Fig. [¥] we plotted the results for the uniform sus- 
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FIG. 5: Uniform spin susceptibility at the wave vector (vr, vr) 
The QMC results are obtained from Ref.— . 



ceptibility Xm=o(0, 0) by using both the bare and full DF 
vertex. The lattice QMC result^^ is shown for compar- 
ison. The calculation is done for U/t = 4.0 and several 
values of temperature. The momentum sum is approxi- 
mated over 32 x 32 points here. Both of these calcula- 
tions reproduce the well known Curie- Weiss law behavior. 
Surprisingly enough, the results for the bare vertex fit the 
QMC results better than that for the momentum depen- 
dent vertex. We believe that this is the finite size effect 
of QMG2^. a. Moreo showed that x becomes smaller 
when increasing the cluster size N. The 4x4 cluster 
calculation result at the same temperature located above 
of that from 8x8 cluster calculation. Therefore the re- 
sults obtained from the full vertex is expected to be more 
reliable. 

The importance of the momentum dependence of the 
DF vertex is more clearly observed in the calculation 
of Xm(7r, tt), see Fig. [S] Again, in this diagram QMC 
results^l are shown for comparison. The same parame- 
ters are used as in Fig. 21 The result from the DF with 
bare vertex does not produce the same results compared 
with QMC solution. Evenmore interesting, with decreas- 
ing temperature the deviation becomes larger. On the 
other hand, the momentum dependent vertex in the DF 
method gives a satisfactory answer. This shows the im- 
portance of the momentum dependence in the DF vertex 
function. Fig. [6] shows the evolution of x against q for 
fixed transfer frequency ujm = 0. The path in momen- 
tum space is shown in the inset. From this diagram we 
can see that x(q, 0) reaches its maximum value at wave 
vector (tt, tt). 

The comparison between the DF and QMC results 
shows the good performance of DF method. The DF 
calculation started from a single site DMFT calculation 
and by introducing an auxiliary field, the non-local in- 
formation is introduced and nicely reproduces the QMC 
results. Our calculation could be done within four hours 
for each value of the temperature on average. In this 
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FIG. 6: xil, 0) vs g at /3i = 2.0, U/t = 4.0 for various q which 
is along the trajectory shown in the inset. 



sense, this method is cheap and reHable compared with 
the more computationally intensive lattice QMC calcu- 
lation. 

Similar as the DF method. Dynamical Vertex Approx- 
imation (DrA)i^ is also based on the two particle local 
vertex. It deals with the lattice fermion directly, with- 
out introducing any auxiliary field. The perturbative na- 
ture of this method ensures its validity at weak-coupling 
regime. Unlike in the DF method, DFA takes the irre- 
ducible two particle local vertex as building blocks. 



7c(l)('''^'''^) =7c(s),,;r('^'^';'^) -Xo{'^;^)Su,^' (24a) 
r-^l^{iy,iy';Q) ^ -f;^l^^„{iy,iy';to) - xo{iy;Q)Su^u' (24b) 



with the spin and charge vertex defined as 7c(s) = 7^^ i 
7^-''. The bare susceptibility is defined as 



Xo{v;uj) = -TGioc{v)Gioc{v + uj) 



(25a) 



T 



Xo[ 



(^., Q) = - - ^ G"(fc, v)G'\k + q,v + uj) (25b) 



And the self-energy is calculated through the standard 
Schwinger-Dyson equation 

^(^) = -^§2 E r/(^' k'- Q)G°{k')G°{k'+Q)G'>{k+Q) 

k',Q 

(26) 
Here, the full vertex F/(fc, k'; Q) is obtained by summing 
all the channel dependent vertices and subtracting the 
double counted diagrams. 

Tf{k,k'-Q) = \l[iV,{v,v';Q)-Vs{v,v'-Q)] 

~[T,{v,u'-uj)^Ts{v,v'-Lo)]\ (27) 



The one particle propagator is given by the DMFT lat- 
tice Green's function where the self energy is purely local 
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FIG. 7: Comparison with the DFA susceptibilities x(0, 0) 
which obtained from both the DMFT lattice Green's func- 
tion (DFA (G")) and the full Green's function (DFA (G)), 
see context for more details. 



G^{k) = l/[ii/ — e{k) — ^{i')], the local Green's function is 
Gioc{v) — l/[iv — A{v) — Yi{v)]. Then the Dyson equation 
gives the lattice Green's function from the self-energy 
fimction G-i = G^^ - S. 

Before presenting the comparison, we take a deeper 
look at the analysis of Eq. (I24p , 

[Xo{i^-.Q)-Xo{v.^)]5..u' (28) 

The second term in the brackets on RHS removes the 
local term from the bare susceptibility. The whole term 
in the brackets then represents only the non-local bare 
susceptibility. In order to compare with the DF method, 
we take the inverse form of Eq. (fT9|) 



r<i,L(^.^';Q) == 7c(l)(^'''''^) 



T 



Y,G\k)G\k + q) (29) 



N 



The above two equations are same except for the last 
term. Since the local DF Green's function Gf^^ is zero, 
the bare DF susceptibility is purely non-local which coin- 
cides with the analysis of DFA Bcthc-Salpeter equation. 
Therefore, it is not surprising that these two methods 
generate similar results. It is not easy to perform a term 
to term comparison between the DF method and DFA 
although the bare susceptibilities have no local term in 
both of these method. The one particle Green's functions 
have different meaning in these two methods. 

The lattice susceptibility within the DFA method is 
obtained by attaching four Green's functions on the ver- 
tex obtained in Eq. (|27p . There are two possible choices 
of the lattice Green's function, one is the DMFT lattice 
Green's function G^, the other one is the Green's func- 
tion G constructed by the non-local self-energy from the 
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FIG. 8: DFA susceptibilities x(7r, tv) at U/t = 4.0. The sus- 
ceptibility are determined from both of the DMFT and full 
lattice Green's function together with the vertex obtained 
from Eq. ^ 



Dyson equation. In Fig. [7]and[8l we presented the DFA 
lattice susceptibility calculated from both the DMFT lat- 
tice Green's function labeled as DrA(G'°) and the full 
Green's function labeled as DFA(G). The DF result from 
the calculation with the full DF vertex is re-plotted for 
comparison. In Fig. [71 the DFA susceptibility calculated 
from the DMFT Green's function (DFA(G°)) is basically 
the same as the DF susceptibility only with some small 
deviation. The results for T/t > 1.0 which are not shown 
here which nicely repeat the DF and QMC results, the 
deviation between the DFA and the DF method becomes 
smaller with the increasing of temperature. The DFA 
susceptibility is calculated from the full Green's func- 
tion (DFA(G')) shows a different behavior at low tem- 
perature regime which reached its maximum value at 
T/t « 0.36. As we know, the Hubbard Model at half fill- 
ing with strong coupling maps to the Heisenberg model, 
X reasches a maximum at T « J where J is the effective 
spin coupling constant given as At'^/U. The calculation 
uses the parameter U/t = 4.0 which is in the intermedi- 
ate coupling regime. Therefore we further calculated the 
lattice susceptibility at U/t = 10.0 which are shown in 
Fig. 1 

When the temperature is greater than 0.4, the DF 
method and DFA (DFA(G°)) generate the similar re- 
sults to the QMC calculation. Reducing the tempera- 
ture further, the QMC susceptibility greatly drops and 
shows a peak around 0.4 which coincides with the be- 
havior of the Heisenberg model. The DF femion and 
DFA susceptibility continuously grows up with the de- 
creasing of temperature. Although the DFA with the 
full Green's function (DFA(G)) shows a peak, it locates 
at T/t — 0.6667 which is larger than the peak position 
of the QMC. And DFA(G) generated a large deviation 
from that of QMC. In this diagram, we only show the 
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FIG. 9: The comparison of the DF resulsts and that of 
QMC for the uniform susceptibility at U/t =10. 4x4 QMC 
results^ also shows the errorbars. 



results of the DF approach for T/t > 0.3 and the DFA 
results for T/t > 0.4. The Bethe-salpeter equation of 
the DFA have a eigenvalue approaching one when fur- 
ther lowering the temperature, which makes the access 
of lower temperature region impossible. 

Fig. [5] shows the results of DFA susceptibilities at 
wave vector (tt, tt). In contrast to the comparison for 
x(0, 0) results, the DFA susceptibility calculated from 
the full Green's function DFA (G) yields better results 
than that from the calculation with the DMFT Green's 
function DFA (G'^). DFA (G) results are almost on top of 
the DF results, the results with DMFT Green's function 
DFA (G°) is large than the DF results. The deviation 
becomes larger at lower temperature. Summarizing, the 
DFA calculation using the full Green's function gener- 
ated the same result as the DF method for xi'^^ ■"") while 
failed to produce x(0, 0) correctly. In contrast, the cal- 
culation with the DMFT Green's function in DFA nicely 
produced the results calculated with the DF method for 
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FIG. 10: The evolution of maximum eigenvalue in spin chan- 
nel against temperature for DF method and DFA. 
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x(0, 0) while generated larger devivation for xl^r, tt) at 
lower temperature regime compared to that from the DF 
method. Together with Fig. [4] and [3 we can see that the 
DF fermion calculation with the full DF vertex generated 
basically the same results for both x(0, 0) and xi'^i ^) 
compared to the results of QMC. 

In both the DF method and the DFA, the operation of 
inverting large matrices is required for solving the Bethe- 
Salpeter equation. Fig. [10] shows the leading eigenvalue 
of Eqns. (|19|1 and ll24|) . As expected, the leading eigen- 
value approaches one with decreasing temperature which 
directly indicates the magnetic instability of 2D system. 
The eigenvalues corresponding to the DF fermion method 
always lies below of that from DFA indicating the better 
convergence of the DF method. When the leading eigen- 
values are closed to one, the matrix inversion in Eqns. 
(fT9|) and ([24]) are ill defined, which prevents the investi- 
gation at very low temperature. 

Concerning the performance of the DF method, we 
also calculated the uniform susceptibility at away half- 
filling. In the strong-coupling limit, the Hubbard model 
is equivalent to the Heisenberg model with coupling con- 
stant J = At'^/U. The consequence of doping is to ef- 
fectively decrease the coupling J, which yields the in- 
creasing behavior of x with doping. The finite size QMC 
calulation^!^ observed a slightly increasing x with very 
small doping at strong interaction or in the low tempera- 
ture region. Here, we did a similar calculation at (3t = 2.5 
and U/t = 4, 10. Since the DF method and the DFA do 
not suffer from the finite size problem. We would ex- 



pect to observe results similar to those of QM C^^i^^ . In 
DFA the suseceptibility is calculated from the DMFT 
Green's function G° and the vertex obtained from Eq. 
(P7| . As shown in Fig. [TT]at U/t = 4.0, the susceptibil- 
ity X slightly increases in the weak dopping region where 
6 is around 0.05, DF fermion results clearly showed such 
behavior, DFA also gave a signal of it. Further doping 
the system, both the DFA and the DF method reproduce 
the decrease with doping as already seen in the QMC. 
With the increasing of interaction, we would expect to 
see the enhancement of this effect, however our calcu- 
lation indicates that such increasing-decreasing behaviro 
dissappear. Both the DFA and the DF method give the 
same decreasing curve which contradict to QMC resul"b2&. 
The results will most likely be further improved by in- 
cluding the higher order vertex or calculating the cluster 
DMFT plus DF/DFASa. 



VI. CONCLUSION 

In this paper, we extended both the DF method and 
DFA to calculate the lattice susceptibility. Both of these 
methods gave equally good results compared with QMC 
calculation at U/t = 4.0. Although they are supposed 
to be weak-coupling methods, at U/t — 10.0 these two 
methods generated right results at high temperature re- 
gion. While both of them failed to reproduce the Heisen- 
berg physics at low temperature. The investigation of the 
lattice susceptibility suffers from hard determined matrix 
inversion problem at low temperature regime. The DF 
methods always generates smaller eigenvalues compared 
to DFA indicating the better convergence. The imple- 
mentation of DF method in momentum space greatly 
improves the calculational speed and makes it easier to 
deal with larger size lattice. 
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